US007124038B2 


d2) United States Patent 

van Dam et al. 


(io) Patent No.: US 7,124,038 B2 

( 45 ) Date of Patent: Oct. 17, 2006 


(54) METHOD AND APPARATUS FOR 

AUTOMATICALLY GENERATING AIRFOIL 
PERFORMANCE TABLES 

(75) Inventors: Cornells P. van Dam, Davis, CA (US); 

Edward A. Mayda, Davis, CA (US); 
Roger Clayton Strawn, Mountain 
View, CA (US) 

(73) Assignee: The Regents of the University of 
California, Oakland, CA (US) 

( * ) Notice: Subject to any disclaimer, the term of this 

patent is extended or adjusted under 35 
U.S.C. 154(b) by 46 days. 

(21) Appl. No.: 11/004,716 

(22) Filed: Dec. 2, 2004 

(65) Prior Publication Data 

US 2005/02461 10 Al Nov. 3, 2005 

Related U.S. Application Data 

(60) Provisional application No. 60/526,686, filed on Dec. 
2, 2003. 

(51) Int. Cl. 

GO IF 19/00 (2006.01) 

(52) U.S. Cl 702/50; 703/1 


(58) Field of Classification Search 702/50 

See application file for complete search history. 

(56) References Cited 

U.S. PATENT DOCUMENTS 


5,544,524 A * 8/1996 Huyer et al 73/147 

6,516,292 Bl* 2/2003 Yahalom 703/9 


* cited by examiner 

Primary Examiner — Michael Nghiem 

Assistant Examiner — Stephen J. Cherry 

(74) Attorney, Agent, or Firm — Park, Vaughan & Fleming 

LLP; Edward J. Grundler 

(57) ABSTRACT 

One embodiment of the present invention provides a system 
that facilitates automatically generating a performance table 
for an object, wherein the object is subject to fluid flow. The 
system operates by first receiving a description of the object 
and testing parameters for the object. The system executes a 
flow solver using the testing parameters and the description 
of the object to produce an output. Next, the system deter- 
mines if the output of the flow solver indicates negative 
density or pressure. If not, the system analyzes the output to 
determine if the output is converging. If converging, the 
system writes the output to the performance table for the 
object. 
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METHOD AND APPARATUS FOR 
AUTOMATICALLY GENERATING AIRFOIL 
PERFORMANCE TABLES 

RELATED APPLICATION 5 

This application hereby claims priority under 35 U.S.C. 
§119 to U.S. Provisional Patent Application No. 60/526,686, 
filed on 2 Dec. 2003, entitled “A General Computational 
Fluid Dynamics-Based Methodology to Automatically Gen- to 
erate Parameterized Airfoil Performance Data,” by inventors 
C. P. van Dam and E. A. Mayda, which is incorporated 
herein by reference 

GOVERNMENT LICENSE RIGHTS 15 

This invention was made with United States Government 
support under NAS A Ames Cooperative Agreement Number 
NCC 2-5485 awarded by Ames Research Center Grant 
Office, Mail Stop 241-1, Moffett Field, Calif. 94035-1000. 20 
The United States Government has certain rights in the 
invention. 

BACKGROUND 

25 

1. Field of the Invention 

The present invention relates to airfoil performance. More 
specifically, the present invention relates to a method and an 
apparatus for automatically generating airfoil performance 
tables. 30 

2. Related Art 

In order to analyze the performance of an airfoil, airfoil 
analysis codes typically make extensive use of two-dimen- 
sional performance data. This two-dimensional airfoil per- 
formance data is typically stored in “C81 airfoil perfor- 35 
mance tables.” C81 tables are text files, which list the lift 
coefficient, drag coefficient and pitching moment coefficient 
as functions of angle of attack for a range of Mach numbers. 
Typically, the angle of attack spans angles from -180° to 
180°, while the Mach number ranges from 0.0 to 1.0. In the 40 
past, tables were created from experimental data and 
researchers’ best estimates as to how the airfoil should 
perform at certain flow conditions. 

Due to the high cost of wind tunnel tests, only a handful 
of airfoil sections have been tested over the entire 360° angle 45 
of attack and Mach number ranges. These data are often 
faired into tables being created for new geometries as a way 
of filling in the gaps where experimental data is not avail- 
able. Due to this process, it is difficult to ensure any level of 
consistency in existing tabular data as far as compressibility 50 
effects, viscous effects or even sectional geometries are 
concerned. 

C81 tables are thus often incomplete, and often contain 
inconsistencies in reference Reynolds numbers, airfoil 
geometry, and the like. These problems make relying on C8 1 55 
tables for design tradeoff studies and geometry modifica- 
tions problematic at best. 

Hence, what is needed is a method and an apparatus for 
automatically generating airfoil performance tables without 
the problems described above. 60 

SUMMARY 

One embodiment of the present invention provides a 
system that facilitates automatically generating a perfor- 65 
mance table for an object, wherein the object is subject to 
fluid flow. The system operates by first receiving a descrip- 


2 

tion of the object and testing parameters for the object. The 
system executes a flow solver using the testing parameters 
and the description of the object to produce an output. Next, 
the system determines if the output of the flow solver 
indicates negative density or pressure. If not, the system 
analyzes the output to determine if the output is converging. 
If converged, the system writes the output to the perfor- 
mance table for the object. 

In a variation of this embodiment, if the output of the flow 
solver indicates negative density or pressure, the system 
adjusts a testing parameter and repeats executing the flow 
solver to produce a new output. The system then determines 
if the new output indicates negative density or pressure. 

In a further variation, if the system mode is “steady-state,” 
adjusting the testing parameter involves decreasing a time 
step for the flow solver. 

In a further variation, if the system mode is “time- 
accurate,” adjusting the testing parameter involves modu- 
lating a Courant-Friedrichs-Levy number limit for the flow 
solver. 

In a further variation, if the system mode is time-accurate, 
adjusting the testing parameter involves adjusting a physical 
time step for the flow solver. 

In a further variation, if the system mode is time-accurate, 
adjusting the testing parameter involves adjusting a pseudo- 
time step for the flow solver. 

In a further variation, if the system mode is steady-state, 
analyzing the output involves analyzing a residual error 
using a regression analysis. 

In a further variation, if the history of the residual error 
indicates numerical divergence, the system sets the system 
mode to time-accurate. 

In a further variation, if the system mode is steady-state, 
analyzing the output involves analyzing a performance 
parameter for the object using a moving average technique, 
wherein the performance parameter includes, but is not 
limited to. lift, drag, and pitching moment. 

In a further variation, if the system mode is time-accurate, 
analyzing the output involves performing a time-correlation 
analysis of a performance parameter history. 

In a further variation, being subject to fluid flow includes 
fluid flow through and/or over the object. 

BRIEF DESCRIPTION OF THE FIGURES 

FIG. 1 illustrates an airfoil performance processing sys- 
tem in accordance with an embodiment of the present 
invention. 

FIG. 2A presents a flowchart illustrating a steady-state 
calculation block in accordance with an embodiment of the 
present invention. 

FIG. 2B presents a flowchart illustrating a time-accurate 
calculation block in accordance with an embodiment of the 
present invention. 

FIG. 3 presents a flowchart illustrating the process of 
varying input parameters in steady-state mode in accordance 
with an embodiment of the present invention. 

FIG. 4 presents a flowchart illustrating the process of 
varying input parameters in time-accurate mode in accor- 
dance with an embodiment of the present invention. 

FIG. 5 presents a flowchart illustrating the process of 
determining if the calculations are complete in accordance 
with an embodiment of the present invention. 

FIG. 6 illustrates curvilinear coordinate transformation 
form the physical domain to the computational domain in 
accordance with an embodiment of the present invention. 
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DETAILED DESCRIPTION 


3 

FIG. 7 illustrates a moving average analysis of a generic 
force or moment coefficient history in accordance with an 
embodiment of the present invention. The following description is presented to enable any 

FIG. 8 illustrates overlapping numerical and physical person skilled in the art to make and use the invention, and 


domains of influence in accordance with an embodiment of 
the present invention. 

FIG. 9 illustrates a physical timestep and CFL max modu- 
lation history in accordance with an embodiment of the 
present invention. 

FIG. 10 illustrates correlation analysis of the lift coeffi- 
cient in accordance with an embodiment of the present 
invention. 

SYMBOLS AND ABBREVIATIONS 

p density 
(j. viscosity 
e energy 

y ratio of specific heats 
x shear stress 
p pressure 

c airfoil chord length 
1 characteristic reference length 
t time 

x,y Cartesian coordinates 
|,r) computational domain spatial coordinates 
u,v velocity components in the x and y directions, respec- 
tively 

a speed of sound, a 2 =yp/ p 
Re Reynolds number, Re=pu 00 c/p 
Pr Prandtl number 
M Mach number, M=\/u 2 +v 2 /a 
a angle of attack 

Cj, C d , C m coefficients of lift, drag and pitching moment 
(about 0.25c), respectively 
dt timestep 

k scaling constant or dissipation coefficient 
e tolerance value 

CFL Courant-Friedrichs-Levy number 
A sample interval used for correlation stencil 
B sample interval being compared to correlation stencil 
n number of values in samples A and B 
j summation index 

s statistical variance variable for correlation analysis 

r correlation coefficient 

R modified correlation coefficient 

SF correlation analysis scaling factor 

<ft,n weight factors used to calculate the correlation factor 

CF correlation factor 

SS steady-state 

TA time-accurate 

MWT maximum wall time allotted per case 
Subscripts 

pseudo denotes pseudo timestep 
ma value pertaining to the moving-average analysis 
max maximum value or upper limit for a variable 
floor absolute lowest value a variable can be assigned 
sf value pertaining to the scaling factor of the correlation 
analysis 

A, B correlation analysis stencil and subinterval 

f farfield values 

oo freestream values 

user user specifiable variable 

x,y Cartesian coordinate directions 

E;,r| curvilinear coordinate directions 


5 is provided in the context of a particular application and its 
requirements. Various modifications to the disclosed 
embodiments will be readily apparent to those skilled in the 
art, and the general principles defined herein may be applied 
to other embodiments and applications without departing 
to from the spirit and scope of the present invention. Thus, the 
present invention is not intended to be limited to the embodi- 
ments shown, but is to be accorded the widest scope 
consistent with the principles and features disclosed herein. 

The data structures and code described in this detailed 
15 description are typically stored on a computer readable 
storage medium, which may be any device or medium that 
can store code and/or data for use by a computer system. 
This includes, but is not limited to, magnetic and optical 
storage devices such as disk drives, magnetic tape, CDs 
20 (compact discs) and DVDs (digital versatile discs or digital 
video discs), and computer instruction signals embodied in 
a transmission medium (with or without a carrier wave upon 
which the signals are modulated). For example, the trans- 
mission medium may include a communications network, 
25 such as the Internet. 

Processing System 

FIG. 1 illustrates an airfoil performance processing sys- 
tem in accordance with an embodiment of the present 
invention. The airfoil processing system includes processor 
106. Processor 106 can generally include any type of com- 
puter system, including, but not limited to, a computer 
system based on a microprocessor, a mainframe computer, a 
digital signal processor, a portable computing device, a 
35 personal organizer, a device controller, and a computational 
engine within an appliance. 

Processor 106 receives airfoil description 102 and pre- 
input file 104 and generates C81 output file 112. Pre-input 
104 file contains the test parameters for testing the airfoil. 
Processor 1 06 can perform two types of calculations: steady- 
state calculations 108, and time-accurate calculations 110. 
Steady-state calculations are typically applicable to conven- 
tional flow conditions such as small angles of attack with 
fully attached flow for which the solution is time-indepen- 
dent. Conversely, time-accurate calculations model the tem- 
poral behavior of time-varying flows. The automation tech- 
niques developed in this specification use both calculation 
techniques either independently or in combination with one 
another. 

50 Airfoil Performance Calculations 

FIG. 2A presents a flowchart illustrating operation of a 
steady-state calculation block 201 in accordance with an 
embodiment of the present invention. The system starts by 
reading a pre-input (PIN) file (step 202). Next, the system 
55 checks for completion of steady-state calculations (step 
204). This check is described in more detail in conjunction 
with FIG. 5 below. If the steady-state calculation is not 
complete, the system generates a flow solver input file (step 
206). Next, the system executes the flow solver using the 
60 input file (step 208). Note that many different flow solvers 
can be used in tins step. Details of a specific flow solver are 
presented below. 

The output of flow-solver 208 is checked for negative 
density and/or pressure (step 210). Note that a negative 
65 density or a negative pressure is not possible in the physical 
world. If the pressure or density is negative, the system 
decreases the pseudo-time step (step 212) and returns to step 
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206 to generate a new flow solver input file. The process of 
decreasing the pseudo-time step is described in more detail 
below in conjunction with FIG. 4. 

If the pressure and density are not negative at step 210, the 
system concatenates the output of the flow solver to an 
output file (step 214). Next, the system analyzes the residual 
history of the force and pitching moments (step 216) and 
increments the run count (step 218). Details of analyzing the 
residual history are presented in more detail below. The 
process then returns to step 204 to determine if the steady- 
state calculations are complete. If so, the system determines 
if satisfactory results were obtained during the steady-state 
calculations (step 220 in FIG. 2B). If not, the system 
switches to time-accurate mode. 

FIG. 2B presents a flowchart illustrating operation of a 
time-accurate calculation block 221 in accordance with an 
embodiment of the present invention. The system first 
determines if the time-accurate calculations are complete 
(step 222). If not, the system modulates the physical time 
step or the Courant-Friedrichs-Levy number upper limit 
(CFL max ) (step 224). Next, the system generates a flow 
solver input file using the modulated physical time step or 
the modulated CFL max (step 226) and executes the flow 
solver using the input file (step 228). Note that the flow 
solver is typically the same flow solver used at step 208. 

The output of flow-solver 228 is then checked for nega- 
tive density and/or pressure (step 210). If the pressure or 
density is negative, the system modulates the physical time 
step or CFL max (step 232) and returns to step 226 to generate 
a new flow solver input file. The process of modulating the 
physical time step or CFL ma;c is presented in detail below in 
conjunction with FIG. 5. 

If the pressure and density are not negative at step 230, the 
system concatenates the output of the flow solver to an 
output file (step 234). Next, the system analyzes the force/ 
moment histories (step 236) and increments the run count 
(step 238). Details of analyzing the force/moment histories 
are presented in detail below. The process then returns to 
step 222 to determine if the steady-state calculation is 
complete. 

If satisfactory steady-state results are obtained at step 220 
or if the time-accurate calculations are complete at step 222, 
the system writes the output file 224 (step 240). Finally, the 
system performs a file cleanup by deleting temporary files 
(step 242). 

Varying Steady-State Input Parameters 

FIG. 3 presents a flowchart illustrating the process of 
varying input parameters in steady-state mode in accordance 
with an embodiment of the present invention. The system 
starts by determining if the pressure or density is negative 
(step 302). If so, the system sets the time step to 0.75 times 
the current time step (step 312). If the pressure and density 
are not negative at step 302, the system determines if the 
current time step is less than the desired time step (step 304). 
If so, the system set the time step to 1 .25 times the current 
time step (step 310) Note, that different multipliers can be 
used in any of these calculations. 

Varying Time-Accurate Input Parameters 

FIG. 4 presents a flowchart illustrating the process of 
varying input parameters in time-accurate mode in accor- 
dance with an embodiment of the present invention. The 
system starts by determining if the pressure or density is 
negative (step 402). If so, the system determines if the 
current CFF max is greater than the CFL max floor (step 406). 
If so, the system sets C 'll max to 0.8 times the current CFL max 
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(step 414). Otherwise, the system sets the time step to 0.75 
times the current time step (step 412). 

If the pressure and density are not negative at step 402, the 
system determines if the current time step is less than the 
5 desired time step (step 404). If so, the system sets the time 
step to 1.333 times the current time step (step 408). Other- 
wise, the system sets CFL max to 1.25 times the current 
CFL (step 410). Note, that different multiplicative param- 
eters can be used in any of these calculations. 

Determining if Calculations are Complete 

FIG. 5 presents a flowchart illustrating the process of 
determining if the calculations are complete in accordance 
15 with an embodiment of the present invention. Note that this 
process relates to both step 204 and 222 above. The system 
starts by determining if the maximum run time or step count 
has been exceeded (step 502). If not, the system then 
determines if the maximum wall time lias been exceeded 
20 (step 504). 

If the wall time lias not been exceeded, the system 
determines if the steady-state has converged or it the time- 
accurate mean-steady-state has been reached depending 
upon the system mode (step 506). If not, the system then 
25 determines if an interrupt or terminate signal has been 
caught (step 508). 

If the system has not caught an interrupt or terminate 
signal, the process returns “no” (step 510). Otherwise, if the 
30 response to any of the above tests has been yes, the system 
returns “yes” (step 512). 

Flow Solver 

The flow solver used in one embodiment of this invention 
35 is ARC2D, a finite-difference two-dimensional compressible 
Reynolds-averaged, Navier-Stokes (RaNS) solver produced 
at NASA Ames Research Center by Pulliam and Steger. The 
ARC2D solver contains many numerical schemes and is 
greatly optimized for computational efficiency, which makes 
40 it an ideal candidate for large parameter studies such as C81 
table generation. This solver is capable of analyzing any 
geometry that can be modeled with a single structured C- or 
O-mesh. Computations can be conducted in either steady- 
state or time-accurate modes. Steady-state calculations uti- 
45 lize space-varying time steps for improved convergence 
rates, and the user is given the power to define mesh and 
time step sequencing scenarios to further enhance compu- 
tational efficiency. Time-accurate solutions can be run with 
a second-order accurate subiteration time-advancing 
50 scheme, in which the user can specify values for parameters 
such as the number of subiterations, physical time step and 
pseudo time step. ARC2D contains an assortment of turbu- 
lence models: Baldwin-Lomax, Baldwin-Barth, Spalart-Al- 
lmaras (several variations), and a two-equation model. 
55 Unfortunately, not all models are compatible with all grid 
geometries and time-accurate calculations. ARC2D also has 
the option to use low-Mach number preconditioning, winch 
is meant to improve the solver accuracy and efficiency as the 
flow approaches incompressible conditions. However, pre- 
go conditioning was not utilized in the current embodiment of 
the present invention. 

In the following sections, ARC2D is discussed in more 
detail with attention paid to the Navier-Stokes equations, 
curvilinear coordinate transformations, various numerical 
65 schemes and formulations, and boundary conditions. The 
one-equation Spalart-Allmaras turbulence model is also 
briefly described. 
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Navier-Stokes Equations 

ARC2D solves the Navier-Stokes equations in strong 
conservation form in an effort to accurately capture shocks. 
The non-dimensionalized equations in Cartesian coordinates 
are shown in Equations. 1, 2 and 3: 


d,Q + d x E+d y F= Re~\d x E y +d y F,) 


(1) 


where 


10 
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pu 


pv 

pu 


pu 2 + p 


puv 
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, E = 


pv 


puv 


pv 2 +p 

. « . 


u(e + p ) . 


v(e + p ) . 


E v = 


■ 0 


■ 0 ■ 

Txx 

, F v = 

Txy 

Txy 


T yy 

. fa 


g4 . 


with 


X» = p(4u x -2iy)/3 

Txy = + K») 


(2) 


15 


20 



Curvilinear Coordinate Transformation 

Body-fitted meshes are required to model arbitrary geom- 
etries. To simplify the differencing schemes, a curvilinear 
coordinate transformation is employed to map the body- 
fitted Cartesian x-y grid points to ^-r| grid points, which lie 
on a uniform mesh with spacings of unit length. FIG. 6 from 
Efficient Solution Methods for The Navier-Stokes Equations 
by T. H. Pulliam shows the transformation of the physical 
domain (shown as a C-grid) to the |-r) computational 
domain. As Equation 8 demonstrates, |,r| and x are functions 
of x, y, and t. 

T-t 


S=i 

•c\=x\(x,y,t) (8) 

Cartesian derivatives can be expressed in the computa- 
tional domain via the chain rule, and the result is given in 
Equation 9: 


d, ' 


1 fa It 

ft 

ft 

= 

0 ft ifr 

ft 



0 ft rj, 

. ft . 


= p(-2u x + 4v y )/3 
fa = uTxx + VT xy + pPr~ l (7 - l)" 1 d x a 2 
gA = UT xy + VT„ + pPF l (y- l)~ l d y a 2 


Finite difference approximations are used to determine the 
transformation metrics (|„ §„ % y , r\„ r\ x , \] y ), which are 
30 defined in Equation 10: 

1 x =Jyr X £,y=-Jx^r-x£ x -y£ y 


The equation of state is used to relate the flow variables 
in Q to the pressure as follows: 


n x =-E%-%=Jx % .-f\,=-x x r\ x -y x r\ y 

35 where J, the metric Jacobian, is 


( 10 ) 


p = (y-l)\e--p(u 2 + v 2 )j 


(4) 


40 


1 

./ = 

W-Vf 


HD 


Generally, variables are non-dimensionalized by their Once the transformation is known, the Navier-Stokes 
freestream quantities as shown in Equation 5. equations can be recast in terms of the generalized curvi- 

linear coordinates as shown in Equation 12: 


P = 



U _ V 

’ a ^ 


e 



45 

(5) 


d r Q + %£+ d„ F = Re- 1 [<%£ v + <9, ft] 


( 12 ) 


If 1 is taken as the reference length (1 is commonly chosen 50 
to be the airfoil chord length), time is non-dimensionalized 
as shown in Equation 6. 


r = 


faco 

T 


(6) 55 


The viscous coefficients are defined in Equation 7 


60 


where 



P 


pU 


pu 


puU +% x p 

q = j~ 1 


, E = J~ l 



pv 


PvU +ftP 


e 


. t/(e+p)-ftp 


F = J ~ 1 


pV 

PUV + T] x p 
PVV + T] yP 


V(e + 


(13) 


P p^la^ (7) 

p - — , Re = 

Poo Poo 


Note that Re is based upon a,, instead of u„, as is 
customary for experimental Reynolds numbers. 


with the contravariant velocities 

U=% r +% x u+% y v and P = Uf +T lx w+1 V v (14) 

The viscous flux tenns from Equation 12 are defined in 
Equation 15: 

£ v =J- 1 g i £ v +| y F v ),f’ v =J- I (»l.A.+U,' F v) 


(15) 
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The t stress terms are also transformed into curvilinear 
coordinates as shown in Equation 16: 

r«=n(4(l^( s + r ua^)-2(| F v s +ri :> ,v„))/3 

V = M<-2 g*«s+rV<^)+4(l > .V|+riyV„))/3 

/ 4 =HT i „+vt i9 ,+/iPr-'(Y-l) _1 (S^ 5 a 2+T t^72a 2 ) 

g^tn^+vt^iPr - 1 (y- 1 )= 1 (l^^+ri^a 2 ) (1 6) 

Numerical Schemes 

ARC2D numerical algorithms are implicit in nature. 
While implicit methods generally require much more com- 
putational work than their explicit counterparts, they benefit 
greatly from superior stability bounds. This advantage is 
particularly advantageous when highly refined grids are 
used, since a conditionally stable explicit scheme may 
require prohibitively small time steps in such cases. 

ARC2D employs a three point implicit time differencing 
scheme of the form shown in Equation 17: 

<?A t d / _n\ A t d - n (17) 

AG =T— F (AQ) + — - F Q + 

1 + at l +ip at 

iT7 A e"" 1 + +-5-^ + ^] 

where 

~n -w+l ~n ~n ~ (18) 

A Q = Q - Q and Q = Q(nAr) 

The 0 and 4> terms can be chosen so that the resulting 
scheme can achieve either first- or second-order accuracy in 
time. The first-order Euler implicit scheme is defined with 
0=1 and 4>=0, and the second-order three-point implicit and 
trapezoidal implicit schemes are attained with 0=1 and §=Vi, 
and 0=1/2 and cj)=0, respectively. For the purposes of this 
study, the Euler implicit scheme was used for steady-state 
calculations while the three-point implicit scheme was used 
for time-accurate dual-timestepping computations. 

When Equation 1 7 is substituted into Equation 12, the flux 
vectors E, E v , F and F v are nonlinear functions of Q. A 
Taylor’s series expansion about Q" is used to linearize the 
equations to second-order accuracy. 

ARC2D approximates the continuous differential opera- 
tors 3 S and 3^ with finite difference operators 8g and 6^ on 
a uniformly spaced computational grid where Ai;=l and 
Ar|=l. If j and k are considered to be grid-point indices, the 
second-order accurate central difference operators are 
defined in Equation 19: 

h u jjr(Uj+ur“j- ij,) 12 and ( 19 ) 

ARC2D also has the capability to estimate derivatives 
with fourth and sixth order accurate central differencing 
schemes. 

When the Navier-Stokes equations are cast in finite dif- 
ference form, the linearized system of equations can be 
arranged into a large banded matrix. Although such a matrix 
is sparse, it is still numerically expensive to solve. To lessen 
this expense, an approximate factorization method is applied 
which transforms the single two-dimensional operator 
matrix into two one-dimensional operators. This factoriza- 
tion is accomplished by neglecting the cross derivative terms 
whose magnitudes scale as the square of the time step. The 
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resulting matrices are in block tridiagonal form and are 
solved much more economically than the unfactored algo- 
rithm. 

While the approximate factorization algorithm is an 
5 improvement over its unfactored counterpart, the implicit 
scheme is still computationally expensive when compared to 
explicit schemes. Further gains in computational efficiency 
can be achieved by reducing the number of variables, or 
blocks, from four (density, x-momentum, y-momentum and 
energy) to three (density, x-momentum and y-momentum), 
which is possible if constant total enthalpy is assumed or if 
similar thermodynamic approximations are used. The 
smaller set of equations can be solved if the thermodynamic 
15 relation replaces the energy equation. 

The blocks in the implicit operators can also be diago- 
nalized making the scheme more efficient. This formulation 
works very well for steady-state calculations since the 
explicit side of the diagonal algorithm is unchanged. There- 
20 fore if the implicit side converges, the diagonal scheme 
yields the same results as would the unfactored scheme but 
at a lower cost. Unfortunately, diagonalization of the 
implicit operator degrades the time accuracy of the scheme 
to first-order in time, which can result in erroneous time- 
25 accurate shock speeds and shock jumps. 

While linear stability analyses show the implicit algo- 
rithm to be unconditionally stable, practical flow phenomena 
such as shocks present stability bounds for the algorithms. 
This problem stems from the fact that numerical schemes 
30 and discrete meshes cannot account for all scales of motion. 
High frequency scales of motion are the most problematic. 
This serious shortcoming can result in inaccurate solutions, 
or worse yet, numerical instability. 

This issue is commonly addressed by adding numerical 
dissipation to the algorithm with a magnitude capable of 
drowning out higher frequencies within a solution without 
being large enough to adversely affect the physical accuracy 
of a calculation. In the early development of implicit finite 
difference codes, constant coefficient dissipation was added 
40 to the algorithm. Later on, schemes based upon upwinded 
differencing came into use, and their one-sided differences 
contain naturally occurring dissipation. Upwinding schemes 
are complex to code and have the disadvantage that the 
amount of dissipation cannot be prescribed explicitly. A 
45 simpler approach is to construct a scheme based upon 
central differences, which allows for a controllable amount 
of artificial dissipation. 

ARC2D employs the nonlinear artificial dissipation 
5Q model developed by Jameson, Schmidt and Turkel. The 
so-called .1ST model combines the second- and fourth-order 
dissipations in one expression shown in Equation 20 (only 
the t, terms are presented here): 


Vf (< r j+U'/j + 1 1 ,t + jt)( e M^QjJc - e !//AfV f AfQjj) 


where the difference operators are 

60 

W QM = 7/2 - qj l / A %qjjc = qj+i,k -qjjc (21) 

= qj,k - qj,k- 1 A n q jz k = qjjt+i - qj,k 

65 

and 
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-continued 

= /t 2 Armax(¥y +u , %■_ u ) 

_ IPj+u - 2 Pj,i +Py-ul 
IPj+l,A + ^Pj,k + Pj-l,k\ 


s'jJ = max(0, k 4 Ai - ejl ) 


The k 2 and k 4 tenns are user specifiable second- and 
fourth-order dissipation coefficients, respectively. The *P 
term is a pressure gradient coefficient that amplifies the 
second order dissipation near shocks. The a JJc terms in 
Equation 20 are spectral radius scaling factors and are 15 
defined as 

o Jtt =\U\+a47W ( 23 ) 

Boundary Conditions 

The numerical formulations of the Navier-Stokes equa- 20 
tions require boundary conditions at body surfaces and 
farfield boundaries. ARC2D contains an assortment of 
boundary condition types including viscous, inviscid, adia- 
batic, and isothermal walls and several characteristic inflow/ 
outflow farfield conditions. For lifting airfoils in subsonic 25 
freestream flow, ARC2D has the capability to apply circu- 
lation-based corrections to the flow variables at the farfield 
boundary. With the freestream velocity components defined 
as u 00 =M co cos(a) and v co =M co sin(a), corrected farfield 
velocities are computed according to Equation 24: 30 


2nr(l - A/£sin 2 ((? - a)) 
pr cos 9 

2nr(\. - Mlvsiir - a )) 


particularly advantageous when unstructured grids are used. 
In contrast, ARC2D and many other structured solvers have 
crude “on-olf ’ logic or short ramps for controlling laminar 
and turbulent boundary layer regions. Simply stated, the 
codes either set the eddy viscosity to zero in the momentum 
equation or zero-out the gradients responsible for the gen- 
eration of eddy viscosity within specified laminar flow 
regions. These regions are usually indicated by grid indices 
along the wall (i.e. a grid point corresponds to a trip 
location). While these implementations ensure the existence 
of the prescribed laminar flow, fully mrbulent boundary 
layer conditions cannot be guaranteed directly downstream 
of a trip. The end-user must exercise care when specifying 
transition locations and verifying that transition does indeed 
occur in the final flow solution. 

Steady-State Run Model 

The steady-state run model runs a single flow calculation 
in the flow solver’s steady-state mode. During the calcula- 
tion, flow solver input parameters are modulated to maintain 
numerical stability, and calculation progress is evaluated 
based upon the findings of a moving-average routine and a 
convergence check of the calculation’s residual history. The 
general behavior of the run model is shown in FIG. 2A 
above. If a steady result is obtained, the calculation will be 
stopped and the results are written to file. However, if a 
divergent residual trend is encountered, the calculation will 
be stopped and a suggestion is made to run the calculation 
in a time-accurate mode. 

For steady-state calculations, numerical stability is most 
affected by the choice of the pseudo-time step. Numerical 
instabilities usually result in the presence of negative density 
or pressure values in the solution. Both are physical impos- 
sibilities, and the flow solver stops the calculation if they 
occur. The automation software detects tins condition and 
reduces the pseudo-timestep. Tins reduction is conducted 
according to Equation 25: 


where the circulation is T= 1 /2M co cC / , the Prandtl-Glauert 
transformation factor is fWl-M„ 2 , and r and 0 are the polar 411 
coordinates for the farfield boundary location relative to the 
airfoil quarter-chord point. The farfield circulation correc- 
tion was only used for steady-state calculations in this study. 

Spalart-Allmaras One-Equation Turbulence Model 45 

The Spalart-Allmaras turbulence model is widely used 
today and was the only model used in this study. This 
one-equation model was assembled using empiricism and 
dimensional analysis, Galilean invariance, and dependence 
on molecular viscosity. It is a largely local model, meaning 50 
that the equation at one point does not depend on the 
solution at other points. However, the model does include a 
wall destruction term for the eddy viscosity that is dependent 
on the distance to the surface. The local nature of the model 
simplifies its implementation on complex grid structures, 55 
including unstructured grids. The model is calibrated in 
two-dimensional mixing layers, wakes, and flat-plate bound- 
ary layers. While these calibrations are well suited to 
attached flow over airfoil geometries, they do not consider 
cases where massively separated flow regions are present, 60 
such as post-stall conditions contained within C 81 tables. 

It should be noted that the treatment of the transitional 
boundary layer region differs between the original model 
and its implementation in ARC 2 D. The Spalart-Allmaras 
model, as first described in the literature, employs a source 65 
term with a small domain of influence to initiate transition, 
or “trip” the flow, in a smooth manner. Tins approach is 


j 4iw _ „ jfCurrent 

ul pseudo ~ K pseudo ' ul pseudo 


where the user defined scaling constant, K. pseudo , can have 
values between 0.0 and 1.0. Its default value is 0.75. After 
every instance of a negative density or pressure condition, 
the pseudo-timestep is reduced. No attempt is made to 
increase the pseudo-timestep. 

Upon the conclusion of each run, which is a small subset 
of the larger calculation, flow solution progress is evaluated 
with a moving-average routine. The averaging routine ana- 
lyzes the lift, drag, and pitching moment coefficient histories 
over the final 50% of the calculation steps or the last 250 
steps whichever is smaller. This parent interval is then 
averaged to determine reference values for the force and 
moment coefficients. The parent interval is broken into 
segments of constant length and averages are computed for 
each segment. This process is repeated for segment lengths 
of 5 to 15 flow solver steps. If all segment averages agree 
with the parent average within a tolerance, the solution is 
deemed steady and complete and the results are written to 
file. A visual explanation of this analysis is offered in FIG. 
7 . 

The tolerances used for the moving average analysis are 
determined as a percentage of the parent interval average 
value. For cases with very small lift, drag or pitching 
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moment results, user definable tolerances are used instead. 
The tolerance is determined from the relationship shown in 
Equation 26. 

e=max (€ usen K__ -parent interval average) (26) 

5 

where is the user defined absolute tolerance and K ma is 
the user defined relative tolerance. 

In addition to checking to see if the force and pitching 
moment histories are behaving in a steady manner, the 
overall trend of the residual history is checked after each 10 
run. This check is accomplished by determining the slope of 
the logarithm of the residual for the most recent run via a 
linear least-squares curve fit. If the slope of this fit is less 
than a user specifiable magnitude (a conservative default 
tolerance is 0.0002), the calculation is deemed convergent; 15 
conversely, if the slope exceeds the tolerance, the calculation 
is considered to be divergent and is stopped, and a sugges- 
tion is made that the calculation be run using the time- 
accurate run model. 

20 

Time-Accurate Rim Model 

The time-accurate run model conducts a single flow 
calculation in a time-accurate dual-timestepping mode. Dur- 
ing the automated calculation, flow solver input variables 
are modulated to maintain numerical stability, and calcula- 25 
tion progress is determined via a correlation analysis of the 
force and moment histories. Steady behavior is also possible 
in time-accurate calculations so the moving average routine 
described above in reference to the steady-state run model is 
conducted in addition to the correlation analysis. The gen- 30 
eral progression of an automated time-accurate calculation 
through the ran model is outlined in FIG. 2A and is 
described above. If all force and moment coefficients pass 
either the correlation or moving average analyses, the cal- 
culation is deemed complete and will be stopped. Final 35 
results are written to file. If the calculation is not complete, 
it will be restarted until suitable results are obtained or one 
of the additional conditions in FIG. 5 is satisfied. 

Two flow solver input parameters, the CFL max number 
and the physical timestep, are modulated in a coupled 40 
fashion to maintain numerical stability during time-accurate 
calculations. Before proceeding, it will be helpful to 
describe the significance of the CFL„ I(M . number in more 
detail. The Courant-Friedrichs-Levy (CFL) condition states 
“[the] full numerical domain of dependence must contain the 45 
physical domain of dependence.” FIG. 8 from Computa- 
tional Gasdynamics by C. B. Laney shows an instance 
where the CFL condition is violated. This is most relevant 
for cases in which supersonic flow exists in regions of small 
grid cells. By restricting the CFL number, the distance waves 50 
travel during a timestep can be limited by decreasing the 
duration of the timestep so that the physical and numerical 
domains of dependence overlap in accordance with the CFL 
condition. Effectively, time is artificially slowed in dense 
grid areas by lowering CFF max . Subiteration techniques, 55 
such as dual-timestepping, can be used to maintain time- 
accuracy even when CFL max is affecting the pseudo- 
timestep in localized grid regions if sufficient subiteration 
convergence is obtained before the solution is advanced in 
time. The number of subiterations required to achieve sec- 60 
ond-order accuracy in time was not investigated during this 
study. A general rule-of-thumb is that the residual should be 
decreased by roughly two orders of magnitude during the 
subiteration process before the solution is advanced in time. 

Modulation of the CFL max parameter is akin to swinging 65 
a double-edged sword. While it can ensure numerical sta- 
bility if it is sufficiently low, excessive lowering of CFL max 
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can force a physically unsteady flow to a steady solution. 
Decreasing the physical timestep can also help to satisfy the 
CFL condition by reducing the distance a wave travels 
during a timestep. The reduction of CFI , max is more desir- 
able since reducing the physical timestep equates to longer 
calculation time due to prolonged flow development. 
Coupled modulation of CFL max and the physical timestep 
was employed to avoid unrealistically steady solutions while 
maintaining acceptable calculation times. 

The modulation methodology of CFL max and the physical 
timestep is depicted in FIG. 4 and is described above. The 
ultimate goal of this modulation process is to guarantee 
numerical stability in the calculation while developing the 
flowfield in as short a time as possible. This is accomplished 
by maintaining stability within the subiteration process via 
CFL max an d stability within the time-advancing process via 
the physical timestep. In the event of an instability indicated 
by a negative density or pressure condition, CFL m< „ is 
reduced while the physical timestep is held constant. This is 
continued until stability is achieved or until CFL max is 
lowered to a value designated CFL max , floor . If CFL max is at 
CFL max , floor, the physical timestep is reduced while CFL max 
is held constant. Generally, timestep reduction is only nec- 
essary during the beginning of a computation due to tran- 
sient starting behavior or high Mach number flow. In the 
event that the preceding run was successful, an increase in 
the physical timestep is attempted if the current timestep is 
smaller than its target value, or CFL max is increased so that 
its effect on the solution is minimized. The modulation 
behavior history for a case at M cc =0.8 and a=60° is shown 
in FIG. 9. 

The automation software operates on the assumption that 
all time-accurate force and moment histories should exhibit 
either steady or repetitive oscillatory behavior. Both condi- 
tions are handled within the calculation progress analyses of 
the time-accurate run model. First, a moving average routine 
is used to determine if the solution is steady. If the force or 
moment coefficient is indeed steady, a value for the coeffi- 
cient is reported and the next coefficient is analyzed. How- 
ever, if the behavior is not steady, the investigation continues 
with a correlation analysis. 

If a solution is unsteady, it should exhibit some sort of 
repeating oscillatory pattern, and to be complete, the solu- 
tion should achieve a mean steady-state meaning that the 
period averages of an oscillatory time history no longer vary 
with time. A correlation analysis is employed to determine 
the time history periods and their averages so that a decision 
regarding the attainment of a mean steady-state can be made. 
In the correlation routine, an interval of calculation steps at 
the end of the force and moment history called a stencil is 
moved back in time and compared to the history it overlaps. 
This is best described by FIG. 10. A correlation factor 1002 
is calculated at each point in time to quantify how well the 
stencil 1006 correlates with the parameter history 1004 it 
overlaps. When the stencil 1006 closely resembles the 
history 1004 at an earlier time, the correlation factor will be 
near unity. A perfect match will result in a correlation factor 
equal to unity. 

Determination of the correlation factor incorporates the 
calculation of a standard correlation coefficient and a scaling 
factor based upon the average value of all the points within 
the stencil 1006 and the average value of the interval to 
which the stencil 1006 is compared. In the following expla- 
nation, values associated with the stencil 1006 are denoted 
by the letter A and the underlying interval by the letter B. To 
form the correlation coefficient one needs the mean values 
(Equations 27 and 28) and variances (Equations 29 and 30) 
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of both samples being compared. The sample covariance 
(Equation 31) is needed as well. 
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The correlation coefficient, r, is given below (Equation 
32). Since negative correlation is of little importance here, 
the correlation coefficient is modified (Equation 33) so that 
spurious results are not obtained later in the analysis. 


r - SAB 

(32) 

Sa^B 


R = max(0, r) 

(33) 


If the stencil 1004 and subinterval behave in similar 
fashions but have different mean values, the correlation 
coefficient will still be near unity. Cases such as these must 
be filtered out with a scaling factor (Equation 34) that 
accounts for the difference between sample means. The user 
specifiable weighting term, is 4.0 by default. 

SF=\-K s y\B-A\ (34) 

The final analysis result is the correlation factor (Equation 411 
35), a weighted combination of the modified correlation 
coefficient and the scaling factor. The user specifiable 
weighting factors <E> and Ilhave default values of 5.0 and 2.0 
respectively. 

45 

CF=SF*"R n (35) 

A value for the correlation factor is calculated each time 
the stencil is moved one flow solver step back in time. An 
example of a correlation factor distribution 1002 is shown in 
FIG. 10. The beginning of periods can be determined from 50 
the peaks of the correlation factor 1002 distribution. If a 
peak has a magnitude within a certain tolerance of unity, its 
location is tagged as the start of a period. Once a user- 
defined number of periods are detected, their averages and 
a total average over all periods are calculated. If all of the 55 
period averages are within a tolerance of the total average 
the force or moment coefficient has reached a mean steady- 
state and the result reported. The tolerance is a user-defined 
percentage of the total average or an absolute user-defined 
tolerance, whichever is larger. 60 

Mixed-Type Run Models 

There exist a large number of flow conditions near the 
onset of stall where the ultimate flow behavior, steady or 
unsteady, is not known before the start of a calculation. For 65 
this reason, combinations of the steady-state and time- 
accurate run models have been implemented. Two such run 


models are included in the automation software. Both start 
calculations as if the steady-state run model was controlling 
them, but if the residual check fails, the computation is 
switched into time-accurate mode. The two rim models 
differ in how this switchover is handled. The first model 
restarts the calculation time-accurately with the solution 
obtained from the last steady-state run. In contrast, the 
second model starts the computation from scratch in time- 
accurate mode if the previous steady-state run exhibited a 
divergent residual trend. 

Choosing the steady-state or time-accurate run models at 
the start of a calculation is analogous to predicting the 
future. This is actually possible for many cases. For instance, 
a case at M„=0.3 and a=2° is likely to be steady while one 
at M^O.3 and a=30° should be unsteady due to massively 
separated flow. This decision or prediction, however, 
requires that the user have some knowledge regarding 
expected flow characteristics. What if the decision is not 
clear because the flow conditions are in a grey area such as 
the onset of stall? Tins is where the mixed-type run model 
is most useful. If the flow is steady, the calculation will finish 
in a steady-state mode, and since steady-state calculations 
are much more efficient than their time-accurate counter- 
parts, a great time savings is realized. If the flow happens to 
be unsteady, at least an attempt was made to attain an 
efficient steady-state solution while not jeopardizing physi- 
cally correct modeling of flow phenomena. 

The biggest disadvantage of the mixed-type run model is 
that it is not foolproof. There are cases that lie on the cusp 
between steady and unsteady flow that may be improperly 
handled by the run model due to the hard coded nature of the 
residual trend check tolerance. A very conservative solution 
to this problem is to run all calculations with the time- 
accurate run model, but this approach would be computa- 
tionally expensive and unnecessary for many cases. For this 
reason, the mixed-type run models have the ability to force 
computations into time-accurate mode based upon their 
angle of attack. The user can specify this time-accurate angle 
of attack range. 

C81 Table Generation 

Until now, only the run models that run single jobs have 
been discussed, but C81 airfoil performance tables require 
that a large number of flow conditions be analyzed. To do so, 
the input file for the automation software must be param- 
eterized according to the desired Mach number and angle of 
attack conditions. The job is then submitted to a computa- 
tional resource and the resulting data gathered and written to 
a C81 table. 

Input files are parameterized with a simple template 
approach. A template for the automation software input file 
is created in which text strings such as AFS_FSMACH and 
AFS_ALPFIA denoted the fields for the freestream Mach 
number and angle of attack. A Perl script then replaces these 
fields with their proper values, write the result out to a new 
file, and use the file to run the automation software. 

OpenPBS 2.3.16, by Veridian Information Solutions, Inc., 
is used to combine three dual AMD Athlon processor 
computers into a compute cluster. A Perl script is written to 
handle job submission to OpenPBS, and OpenPBS coordi- 
nated the submission of jobs to specific compute nodes. 
Each case is treated as a serial job that is assigned to its own 
CPU. As a result of this submittal approach, the generation 
of a C81 table is highly scalable, and the time required to 
generate the table is directly related to the number of 
processors available for computation. 
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Every supercomputing center is likely to have its own 
batch processing software and rules governing user eti- 
quette. For this reason, no generalized job submittal soft- 
ware presented for the present invention. Hie end user of the 
automation software is responsible for job submittal, but this 5 
task can be accomplished with existing tools that handle 
parameter studies and job submission to many types of 
computational resources. 

With C81 tables normally consisting of 500 or more 
M^-a pairings, tools are needed for determining the status to 
of jobs and harvesting final data into C81 tables. Once again, 
Perl is used to create these tools. This software allows the 
user to query the C81 table cases for completeness. If a job 
needs to be restarted, its flow parameters and file archive 
location can be written to a tab-delimited text file. The user 15 
can view the progression of restart cases listed in these files 
by plotting the force/moment, residual, physical timestep, 
and CFL max histories. If it appears that the calculation has 
progressed to the point of being complete but the automation 
software indicates otherwise, the user can override this by 20 
manually determining results using a correlation analysis 
with relaxed tolerances or a brute force time averaging 
routine. Lastly, the results from all cases are written to a C81 
table file. 

The foregoing descriptions of embodiments of the present 25 
invention have been presented for purposes of illustration 
and description only. They are not intended to be exhaustive 
or to limit the present invention to the forms disclosed. 
Accordingly, many modifications and variations will be 
apparent to practitioners skilled in the art. Additionally, the 30 
above disclosure is not intended to limit the present inven- 
tion. The scope of the present invention is defined by the 
appended claims. 

What is claimed is: 

1. A method for automatically generating a performance 
table for an object, wherein the object is subject to fluid flow, 
comprising: 

receiving a description of the object, wherein the object is 4Q 
a physical object and wherein the description is a 
curvilinear coordinate transformation of the physical 
object; 

receiving testing parameters for the object; 

executing a flow solver using the testing parameters and 45 
the description of the object to produce an output; 

determining if the output of the flow solver indicates 
negative density or pressure; 

if not, 

analyzing the output to determine if the output is 50 
converging, and 

if converged, writing the output to the performance 
table for the object; and 

using the performance table to analyze the performance of 
the object when the object is subjected to the fluid flow. ^ 

2. The method of claim 1, wherein if the output of the flow 
solver indicates negative density or pressure, the method 
further comprises: 

adjusting a testing parameter; 

executing the flow solver again to produce a new output; 
and 

detennining if the new output indicates negative density 
or pressure. 

3. The method of claim 2, wherein if a system mode is 65 
“steady-state,” adjusting the testing parameter involves 
decreasing a time step for the flow solver. 
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4. The method of claim 2, wherein if a system mode is 
“time-accurate,” adjusting the testing parameter involves 
modulating a Courant-Friedrichs-Levy number limit for the 
flow solver. 

5. The method of claim 2, wherein if a system mode is 
time-accurate, adjusting the testing parameter involves 
adjusting a physical time step for the flow solver. 

6 . The method of claim 2, wherein if a system mode is 
time-accurate, adjusting the testing parameter involves 
adjusting a pseudo-time step for the flow solver. 

7. The method of claim 1, wherein if a system mode is 
steady-state, analyzing the output involves analyzing a 
residual error using a regression analysis. 

8 . The method of claim 7, wherein if the history of the 
residual error indicates numerical divergence, the method 
involves setting the system mode to time-accurate. 

9. The method of claim 7, wherein if a system mode is 
steady-state, analyzing the output involves analyzing a per- 
formance parameter for the object using a moving average 
technique, wherein the performance parameter includes one 
of a lift, a drag, and a pitching moment. 

10. The method of claim 1, wherein if a system mode is 
time-accurate, analyzing the output involves performing a 
time-correlation analysis of a performance parameter his- 
tory. 

11. The method of claim 1, wherein being subject to fluid 
flow involves at least one of: 

fluid flow through the object; and 

fluid flow over the object. 

12. A computer-readable storage medium storing instruc- 
tions that when executed by a computer cause the computer 
to perform a method for automatically generating a perfor- 
mance table for an object, wherein the object is subject to 
fluid flow, the method comprising: 

receiving a description of the object, wherein the object is 
a physical object and wherein the description is a 
curvilinear coordinate transformation of the physical 
object; 

receiving testing parameters for the object; 

executing a flow solver using the testing parameters and 
the description of the object to produce an output; 

detennining if the output of the flow solver indicates 
negative density or pressure; 

if not, 

analyzing the output to determine if the output is 
converging, and 

if converged, writing the output to the performance 
table for the object; and 

using the perfonnance table to analyze the performance of 
the object when the object is subjected to the fluid flow. 

13. The computer-readable storage medium of claim 12, 
wherein if the output of the flow solver indicates negative 
density or pressure, the method further comprises: 

adjusting a testing parameter; 

executing the flow solver again to produce a new output; 
and 

detennining if the new output indicates negative density 
or pressure. 

14. The computer-readable storage medium of claim 13, 
wherein if a system mode is “steady-state,” adjusting the 
testing parameter involves decreasing a time step for the 
flow solver. 

15. The computer-readable storage medium of claim 13, 
wherein if a system mode is “time-accurate,” adjusting the 
testing parameter involves modulating a Courant-Friedrichs- 
Levy number limit for the flow solver. 
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16. The computer-readable storage medium of claim 13, 
wherein if a system mode is time-accurate, adjusting the 
testing parameter involves adjusting a physical time step for 
the flow solver. 

17. The computer-readable storage medium of claim 13, 
wherein if a system mode is time-accurate, adjusting the 
testing parameter involves adjusting a pseudo-time step for 
the flow solver. 

18. The computer-readable storage medium of claim 12, 
wherein if a system mode is steady-state, analyzing the 
output involves analyzing a residual error using a regression 
analysis. 

19. The computer-readable storage medium of claim 18, 
wherein if the history of the residual error indicates numeri- 
cal divergence, the method involves setting the system mode 
to time-accurate. 

20. The computer-readable storage medium of claim 18. 
wherein if a system mode is steady-state, analyzing the 
output involves analyzing a performance parameter for the 
object using a moving average technique, wherein the per- 
formance parameter includes one of a lift, a drag, and a 
pitching moment. 

21. The computer-readable storage medium of claim 12, 
wherein if a system mode is time-accurate, analyzing the 
output involves performing a time-correlation analysis of a 
performance parameter history. 

22. The computer-readable storage medium of claim 12, 
wherein being subject to fluid flow involves at least one of: 

fluid flow through the object; and 

fluid flow over the object. 

23. An apparatus for automatically generating a perfor- 
mance table for an object, wherein the object is subject to 
fluid flow, comprising: 

a receiving mechanism configured to receive a description 
of the object, 

wherein the receiving mechanism is further configured to 
receive testing parameters for the object; 

an executing mechanism configured to execute a flow 
solver using the testing parameters and the description 
of the object to produce an output; 

a determining mechanism configured to determine if the 
output of the flow solver indicates negative density or 
pressure; 

a analyzing mechanism configured to analyze the output 
to determine if the output is converging; and 
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a writing mechanism configured to write the output to the 
performance table for the object. 

24. The apparatus of claim 23, further comprising: 

an adjusting mechanism configured to adjust a testing 

5 parameter; 

wherein the executing mechanism is further configured to 
execute the flow solver again to produce a new output; 
and 

wherein the determining mechanism is further configured 

to to determine if the new output indicates negative den- 
sity or pressure. 

25. The apparatus of claim 24, wherein if a system mode 
is “steady-state,” adjusting the testing parameter involves 
decreasing a time step for the flow solver. 

15 26. The apparatus of claim 24, wherein if a system mode 

is “time-accurate,” adjusting the testing parameter involves 
modulating a Courant-Friedrichs-Levy number limit for the 
flow solver. 

27. The apparatus of claim 24, wherein if a system mode 

20 is time-accurate, adjusting the testing parameter involves 

adjusting a physical time step for the flow solver. 

28. The apparatus of claim 24, wherein if a system mode 
is time-accurate, adjusting the testing parameter involves 
adjusting a pseudo-time step for the flow solver. 

25 29. The apparatus of claim 23, wherein if a system mode 

is steady-state, analyzing the output involves analyzing a 
residual error using a regression analysis. 

30. The apparatus of claim 29, wherein if the history of the 
residual error indicates numerical divergence, the method 

30 involves setting the system mode to time-accurate. 

31. The apparatus of claim 29, wherein if a system mode 
is steady-state, analyzing the output involves analyzing a 
performance parameter for the object using a moving aver- 
age technique, wherein the performance parameter includes 

35 one of a lift, a drag, and a pitching moment. 

32. The apparatus of claim 23, wherein if a system mode 
is time-accurate, analyzing the output involves performing a 
time-correlation analysis of a performance parameter his- 
tory. 

40 33. The apparatus of claim 23, wherein being subject to 

fluid flow involves at least one of: 

fluid flow through the object; and fluid flow over the 
object. 



